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1.  Introduction 


ALEGRA  is  a  multipurpose  computational  tool  for  modeling  a  variety  of  mechani¬ 
cal  and  electromagnetic  phenomena  and  their  interplay.  This  multiphysics  capabil¬ 
ity  is  a  key  feature  of  ALEGRA  and  the  result  of  many  years  of  multidisciplinary  ef¬ 
fort.  However,  it  carries  with  it  a  necessarily  high  cost  for  validation  and  verification 
(V&V).  Eortunately,  V&V  efforts  can  take  advantage  of  the  compartmentalization 
ingrained  in  the  architecture  of  ALEGRA.  In  other  words,  various  modules  of  ALE¬ 
GRA  can  be  used  without  others  when  necessary.  Therefore,  the  V&V  procedures 
can  be  compartmentalized  as  well. 

In  this  report,  we  study  the  reliability  of  ALEGRA’s  transient  magnetics  module  for 
the  static  limit,  and  we  explore  the  quasi-static  evolution  through  which  a  steady 
state  solution  is  reached  for  a  particular  problem  of  interest.  We  begin  in  Section  2 
by  reviewing  the  system  of  equations  used  to  describe  quasi-static  magnetization. 
Eor  the  reader’s  convenience,  an  analytic  solution  is  presented  here  in  summary 
form  for  the  problem  of  interest  in  the  project.  In  Section  3,  the  full  derivation  of 
this  solution  is  presented  in  detail.  This  section  may  be  skipped  in  the  first  reading 
of  the  report.  Section  4  provides  some  important  information  on  the  system  of  units 
to  be  used.  Section  5  describes  the  ALEGRA  code.  Section  6  describes  the  config¬ 
uration  of  the  ALEGRA  simulations.  In  Section  7,  the  qualitative  behavior  of  the 
computed  solution  over  time  is  described  and  illustrated  graphically.  The  quantita¬ 
tive  verification  analysis  appears  in  Section  8,  followed  by  conclusions  in  Section 
9. 


1 


2.  The  Quasi-Static  Theory  of  Polarization  of  Metals 


In  the  theoretical  part  of  this  project  we  follow  the  classical  textbook  of  theoretical 
physics.^ 

2.1  Master  System 

The  analysis  of  quasi-statics  is  based  on  the  following  master  system 


IdB^ 
c  dt 


(1) 


(2) 


/  =  a^^Ej 


(3) 


5*  =  {W)  (4) 

These  bulk  partial  differential  equations  (PDEs)  should  be  augmented  with  the 
boundary  conditions 


=  0,  [H,]tB  =  0  (5) 

and  appropriate  conditions  at  infinity  as  well  as  with  the  corresponding  initial  con¬ 
ditions. 

In  the  system  1-5,  z’'  and  t  are  the  spatial  (Eulerian)  coordinates  and  time;  Ei,  Hi 
and  5*  are  the  electric  and  magnetic  field  and  magnetic  induction,  respectively;  J* 
is  the  electric  current  density  of  free  charges,  c  is  the  speed  of  light  in  vacuum,  cr*-^ 
is  electrical  conductivity. 

Subsystem  1,  2  comprises  a  correspondingly  reduced  classical  Maxwell  system  for 
electromagnetic  field.  Equation  1  is  Earaday’s  law.  Equation  2  is  Ampere’s  law, 
where  displacement  currents  have  been  neglected.  Equation  3  represents  Ohm’s  law 
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in  its  the  simplest  linear  version.  The  funetion  5*  (H^)  is  the  so-ealled  eonstitutive 
equation  for  the  magnetie  field.  The  eonstitutive  relations  3  and  4  should  be  derived 
from  experiments.  In  order  to  avoid  major  physieal  errors  these  equations  should  be 
taken  in  sueh  a  way  that  the  total  master  system  1-4  satisfies  the  First  and  Seeond 
laws  of  thermodynamies.  These  demands  entail  quite  strong  a  priori  eonstraints  on 
the  eonstitutive  equations. 

In  the  boundary  eonditions  5,  rii  and  r*  are  the  normal  and  the  tangent  veetors  to 
the  diseontinuity  boundaries. 

2.2  Notation 

The  metries  eo-  and  eontra-variant  tensors  Zij,  of  the  Eulerian  eoordinate  sys¬ 
tem  are  used  for  lowering  and  raising  (“juggling”)  the  indexes,  and  for  defining 
the  eovariant  differentiation  V*  with  respeet  to  the  eoordinates  z"]  z"^^^  is  the  so- 
ealled  eovariant  Levi-Civita  skew-symmetrie  tensor.  Using  tensorial  notation  per¬ 
mits  one  to  present  all  the  equations  in  the  universal  eovariant  (i.e.,  eoordinate- 
independent)  form.  The  same  goal  ean  be  reaehed  with  the  help  of  the  so-ealled 
direet  or  eoordinate-independent  form,  whieh  in  the  pre-Einsteinian  era  was  the 
most  widespread  in  the  eleetromagnetie  literature.  Eurther  details  on  tensorial  alge¬ 
bra  and  ealeulus  ean  be  found  in  a  plethora  of  textbooks,  ineluding  the  most  popular 
and  praetieal.^’^ 

We  assume  the  eoordinate  system  z''  to  be  Cartesian.  In  the  Cartesian  eoordinate 
system,  the  metries  tensors  numerieally  eoincide  with  the  eorresponding  values  of 
the  Kronnekker  symbol  5^^,  ,  5*.  The  Eevi-Civita  tensor  has  the  main  eom- 

ponent  equal  to  1.  Juggling  of  the  indexes  does  not  ehange  the  numerieal  values 
of  tensor  eomponents.  Nonetheless,  we  will  place  all  of  the  indexes  in  the  covariant 
manner,  thus  following  the  demand  of  mathematical  aesthetics  and  simplifying  the 
passage  between  different  coordinate  systems  when  necessary. 

2.3  Vector  Potential 

The  AEEGRA  code  uses  the  vector  potential  Aj.  The  vectors  Aj  and  Hi  are  inter¬ 
connected  by  the  covariant  differential  relation 


(6) 
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With  a  known  spatial  distribution  of  the  vector  potential  potential  {z,  t)  one  can 
immediately  and  explicitly  recover  the  magnetic  field  using  the  formula  6.  The 
opposite  is  not  true.  First  of  all,  to  recover  the  potential  Ak  with  the  help  of  a  known 
field  iF*  one  has  to  solve  the  boundary  value  problem,  based  on  the  PDF  6.  Equa¬ 
tion  6  alone  is  insufficient  for  solving  this  inversion  problem.  In  order  to  make  the 
inversion  unique  one  has  to  add  to  6  one  of  the  calibration  conditions: 


(a)  the  Coulomb  gauge 


VMi  =  0 


(7) 


or 

(b)  the  (relativistically  invariant)  Lorentz  gauge 

-^  +  VM,  =  0  (8) 

c  at 

where  {z,  t)  is  the  scalar  potential  of  the  electric  field  Ek. 

2.4  System  to  be  Solved 

We  proceed  with  the  analysis  of  linear  isotropic  materials  that  are  free  of  bulk  elec¬ 
tric  charge  or  material  magnetization,  for  which 


a 


ij 


az  ■ 


(9) 


and 


=  pH\ 


(10) 


where  /i  is  the  magnetic  permeability. 

Eliminating  J*  between  2  and  3,  we  get 

A^^V,Hk  =  —E^  (11) 

c 


Eliminating  E^  between  1  and  11,  we  get 


4 


(12) 


(95* 


- — 

Ana  '  ' 


or 


Ana^  dW 
dt 


VjV^5* 


(13) 


The  “parabolic”  equation  13  should  be  combined  with  the  boundary  conditions. 
In  the  equilibrium,  time-independent,  configuration  we  get  the  following  system: 


Vjyw  =  0  (14) 

whereas  the  boundary  conditions  can  be  rewritten  as  follows: 


[/i5,]+n'  =  0,  [5,]+5  =  0  (15) 

If  the  magnetic  field  H°  at  infinity  is  given,  the  system  15  has  to  be  analyzed  under 
the  additional  boundary  condition  at  infinity 


5*  =  5*  ^  =  5*°  (16) 

|z|^oo 


2.5  Elliptic  Inclusion 

There  are  few  exact  2-dimensional  (2-D)  and  3-dimensional  (3-D)  solutions  of  the 
system  14-16.  One  of  them  is  a  solution  for  an  elliptic  inclusion  in  an  infinite 
isotropic  matrix,  in  particular,  in  vacuum.  This  solution  is  described  below. 

Consider  an  ellipse  with  the  semi-axes  a  and  b  coinciding  with  the  Cartesian  axes 
and  as  shown  in  Fig.  1.  We  assume  that  the  elliptical  domain  is  filled  with  a 
linear  isotropic  substance  with  magnetic  permeability  fx.  We  then  assume  that  the 
ellipse  is  immersed  in  the  unbounded  space  in  which  there  is  a  uniform  magnetic 
field  5*°. 
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Fig.  1:  Schematic  diagram  (not  to  scale)  showing  problem  to  be  considered. 


If  there  is  an  elliptical  inclusion,  the  otherwise  uniform  field  =  H^°  will  change. 
The  changes  are  particularly  strong  inside  the  ellipse  and  in  its  vicinity.  At  infinity, 
the  newly  generated  field  if*  approaches  its  original  value,  fi*°. 

2.6  Summary  of  Solution 

Lengthy  calculations  show  that  the  field  if*  inside  the  ellipse  is  uniform,  and  it  is 
given  by  the  relationship 


TTl 

inside 
tt2 
inside 

It  should  be  noticed  that  this  is  a  simple,  elegant,  closed-form  expression  that  is 
very  easily  evaluated.  Because  it  implies  that  the  field  inside  the  ellipse  is  uniform, 
it  can  even  be  used  for  quick  visual,  qualitative  inspection  of  computed  solutions. 

The  field  ff*  outside  of  the  ellipse  is  given  by  the  more  complicated  relationship 


a-\-b 
a-\-fj,b 
b-\-a  TT 
LLa-\-b  ‘ 


(17) 
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Hi 


outside 


outside 


Hi 


(18) 


outside 


outside 


1  d^B  (a+^)(/^— 1)  TTo  I  J_  d^B  {a+b){fi—l)  j^o  _i_  uo 

277  dz^dz^  a-^fib  ^  '  277  dz^dz"^  ^a+b  2  1 

1  d^B  (a+^(/^— 1)  TTO  _i_  J_  d^B  (a+fo)(/^— 1)  tto  _i_  tto 

277  dz‘^dz^  a+^b  ^  '  277  dz'^dz'^  /ia+6  2  '  2 


where  0  is  the  logarithmie  potential  of  the  ellipse,  given  by  the  relationship: 


(19) 


ell 


The  full  derivation  of  this  solution  is  deseribed  in  detail  in  the  following  seetion, 
whieh  is  ineluded  only  for  eompleteness.  This  seetion  may  be  passed  over  in  the 
first  reading. 

3.  Full  Derivation:  Magnetization  of  Isotropic  Elliptical  Inclusion 

This  problem  was  analyzed  by  many  outstanding  mathematicians  and  physicists 
working  on  it  since  Newton’s  times.  First,  it  was  focused  on  various  problems  of 
gravitation  and  cosmology.  Later  on,  in  addition  to  mathematicians,  the  beauty  of 
this  problem  hypnotized  researchers  exploring  electromagnetism,  theory  of  elastic¬ 
ity  among  others.  Monographs ^  and  papers®’’  give  a  picture  of  various  applica¬ 
tions  of  this  problem.  In  the  presentation  below  we  use  the  approach  of  Grinfeld,^ 
which  combines  some  techniques  of  Sretensky"^  and  Eshelby.®’’ 

We  will  be  seeking  a  solution  Hi  =  rji^  (z)  inside  the  ellipse  in  the  following  form: 


rii_  (z)  =  Ki 


(20) 


The  solution  in  the  form  20  automatically  satisfies  the  PDF  14. 

For  the  solution  outside  of  the  elliptic  inclusion,  the  solution  has  the  following  form: 


r7,+  (z)  = V,  V,0+ + 


(21) 
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where  0  is  the  logarithmic  potential  of  an  ellipse. 

The  potential  0  (z)  is  defined  in  Equation  19  and  satisfies  the  PDEs 

ViV*0  =  -27r,  2;  e  Ueiu  (22) 

ViV*0  =  0,  e  R^/oZell 

The  solution  in  the  form  21  automatically  satisfies  the  PDE  and  the  condition  at 
infinity  16. 

The  potential  inside  the  ellipse  0_  (z)  is  described  by  the  quadratic  form 

Q-{z)  =  C-^Y,^Rz^  (23) 

where  C  is  a  constant,  and  Yij  is  a  symmetric  tensor  depending  on  the  shape  of  the 
ellipse  and  described  below. 

We  turn  now  to  the  boundary  conditions  15.  Eirst  of  all,  we  get 

V,V,0+  =  V,V,0_  +  [V,V,0]+  (24) 

Using  the  so-called  compatibility  condition 

[VjVj0]i')  =  27rninj  (25) 

and  the  relationship  23,  we  can  rewrite  24  as  follows: 

ViVj0+o  =  -Yij  +  2Tininj  (26) 

Using  21  and  26,  we  can  rewrite  the  boundary  values  ^7i+o(^)  as  follows: 

+  H°  =  (Yi,  -  2Tinin,)  +  (27) 
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The  solution  0+(2:)  outside  of  the  ellipse  is  much  more  complex.  Fortunately,  hav¬ 
ing  known  the  solution  inside  the  inclusion,  that  is,  the  constants  C  and  Yij,  we  can 
uniquely  recover  the  potential  outside  of  the  inclusion.  Using  the  boundary  condi¬ 
tions  18,  19,  we  can  find  all  of  these  constants. 

Again  by  luck,  permitting  one  to  find  the  exact  solution  is  the  fact  that  the  boundary 
conditions  18,  19  require  one  to  know  the  field  0±(2:)  not  everywhere  but  only  at 
the  ellipse’s  interface  0+o. 

We  can  rewrite  the  boundary  conditions  14,  15  as  follows: 


+  H°-  Ki)  Y  =  0  (28) 

and 


-  /i+  [W^  {Yi^  -  27iSij)  +  H°)}Y  =  0  (29) 

Since  conditions  28,  29  should  be  satisfied  for  arbitrary  orientations  of  the  vectors 
Ti  and  n\  the  system  28,  29  implies 


Ki  -  W% 


H° 


(30) 


and 


fx.Ki  -  /r+  {Yij  -  27i6ij)  (31) 


Using  30,  we  can  rewrite  31  as  follows: 


(/i_  —  /i+)  Ki  -|-  27rmi  =  0 


(32) 


as  implied  by  the  chain 
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0 


(33) 


Ki  -  {Yi,  W^'  -  27t6^jA^)  = 

fi.Ki  -  /i+  {Ki  -H°-  27T6ijW^)  = 

Ki  (/i_  —  /i+)  +  2'!iWi  = 


In  the  case  of  a  body  immersed  in  vacuum  (/i+  =  1  for  Gaussian  units)  the  relation¬ 
ship  32  reads 


(/i_  —  1)  Ki  +  2TiWi  =  0  (34) 

Summarizing,  we  have  arrived  at  2  linear  vectorial  Equations  30,  34  with  2  un¬ 
knowns  Ki  and  Wi.  Excluding  Wi  between  these  2  equations,  we  get  a  single  vec¬ 
torial  equation: 


(^4*  +  K'^  =  H°  (35) 

After  solving  Equation  35  with  respect  to  we  can  find  Wi  from  the  equation 

=  (36) 

271 

thus  finishing  the  determination  of  the  exact  solution. 

In  the  coordinate  system  whose  axes  coincide  with  the  axes  of  ellipse,  the  geometric 
tensor  Yki  has  the  following  components 


\\Ym 


Yu 
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b  0 

3^21 

Y22 

a  +  h 

0  a 

Thus,  we  can  rewrite  the  system  37  as  follows 


(37) 
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0 


X 


1  + 


b 


(/i-i) 


0 


Ki 

HI 

to 

H2 

1  + 


a 

(i-\-b 


(/i-i) 


or 


a+/j6 

a-\-b 


0 


0 


b-\-{ia 

a-\-b 


HI 

to 

(38) 


The  system  38  has  the  following  solution: 


K, 

K2 


(i-\-b 

a-\-^b 

b-\-a 

fj,a-\-b 


HI 


(39) 


and 


Wi 

1  (a+b)(/j,—l)  TTo 
‘I'K  a-\-^b  1 

W2 

1  {a+b)(^—l)  tto 

277  fia-\-b  2 

(40) 


thus  reprodueing  Equation  17.  Then,  inserting  the  relationships  40  of  Wi  and  W2 
in  Equation  21,  we  arrive  at  Equation  18. 

4.  The  Master  System  in  the  SI 

So  far  our  analysis  relied  on  the  usage  of  the  Gauss  system  of  eleetromagnetie  units. 
On  the  other  hand,  the  “magnetohydrodynamics”  (MHD)  module  of  the  AEEGRA 
code  uses  only  Systeme  International  (SI)  units.  Each  of  the  systems  has  its  own 
pros  and  contras.  Before  discussing  the  AEEGRA  results,  let  us  first  rewrite  the 
master  system  1-5  using  SI  units.  We  get 


dB^ 


(41) 


=  f 


(42) 
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In  the  case  of  a  linear  isotropic  medium  that  is  free  of  bulk  electric  charge  or  mate¬ 
rial  magnetization,  we  still  have  the  relationships 


=  /rif* 


(43) 


and 


J*  =  aE^ 


(44) 


Boundary  conditions  are  given  by 


=  0,  [i/.]^r*  =  0 


(45) 


We  note  that  in  SI  units,  the  vacuum  permeability  is  given  by  the  fundamental  con¬ 
stant  /To  =  dvr  X  10“^  N/A^.  For  regions  exterior  to  the  ellipse,  /r  =  /io- 

5.  ALEGRA  Formulation  

The  system  41-43  is  a  form  of  the  equations  solved  by  the  “transient  magnetics” 
module  of  the  ALEGRA  code.^’^°  In  this  module,  the  time  evolution  of  electro¬ 
magnetic  fields  in  a  stationary  conducting  medium  is  computed  in  quasi-static  fash¬ 
ion.  The  system  can  be  obtained  from  Maxwell’s  equations  by  use  of  the  so-called 
“magneto-quasi-static”  approximation.  ^  ^ 

Electromagnetic  radiation  and  separation  of  charges  are  not  possible  in  the  medium 
under  this  approximation.  This  can  be  understood  by  the  absence  of  Maxwell’s 
displacement  current  from  Equation  42,  and  the  absence  of  Gauss’s  law  (relating 
the  electric  displacement  to  the  charge  density)  from  the  system  completely.  The 
divergence-free  character  of  the  magnetic  field  is  implied  in  Equation  45. 

The  medium  is  assumed  to  be  both  stationary  and  electrically  neutral  in  the  “tran¬ 
sient  magnetics”  module.  Mechanical  and  electromagnetic  forces  are  not  applied, 
and  the  thermodynamic  state  (internal  energy)  of  the  material  is  not  updated.  An 
implicit  linear  solver  (see  Section  6)  is  used  to  solve  the  system  with  an  unstruc¬ 
tured  2-D  or  3-D  finite-element  discretization,  and  evolve  the  solution  forward  in 
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time.  (Only  quadrilaterals  are  eonsidered  here,  although  triangles  are  possible  for 
the  “transient  magneties”  module.)  Speeialized  “algebraie  multigrid”  methods  are 
used  to  ensure  suffieient  sealability  of  the  solve  operation  for  large  meshes  (more 
than  10^  elements)  and  parallel  proeessing. 

ALEGRA  is  equipped  to  handle  -  and  has  been  used  sueeessfully  for  -  a  mueh 
broader  elass  of  problems.  These  inelude  deforming  solid  and  hydrodynamie 
media,  meehanieal  and  eleetromagnetie  forees,  and  adiabatie  and  Ohmie  sourees 
of  heating,  via  its  “magnetohydrodynamies”  (MHD)  module.  The  MHD  physies 
module  eouples  the  transient  magneties  eapabilities  with  the  eontinuum  meehanies 
eapabilities  that  form  ALEGRA’s  basis,  enhaneed  to  ineorporate  eleetromagnetie 
induetion  and  eonservation  of  magnetie  flux.  Thus,  the  insights  gained  from  this 
study  should  apply  to  those  problems,  sinee  the  transient  magneties  eapabilities  are 
used  in  eoupled  fashion  in  those  types  of  ealeulations. 

To  obtain  a  simple  transient  magneties  simulation  in  AEEGRA,  with  the  limitations 
deseribed  above,  it  is  only  neeessary  to  use  the  “transient  magneties”  keyword  in 
the  input  as  the  name  of  the  physies  module  to  be  employed.  AEEGRA  solves  either 
the  system  41-43,  or  the  system  reeast  in  terms  of  the  veetor  potential  Ai,  depending 
on  the  eonfiguration  of  the  problem.  Eor  the  problem  eonsidered  here,  the  geometry 
is  sueh  that  AEEGRA  solves  the  system  in  terms  of  the  out-of-plane  eomponent  of 
the  veetor  potential,  whieh  we  refer  to  as  A.  The  magnetie  flux  density  Bi  is  then  a 
derived  variable. 

6.  Simulation  Development 

Two-dimensional  AEEGRA  transient  magneties  simulations  were  run  for  the  prob¬ 
lem  deseribed  in  Seetion  2  and  Eig.  1.  In  partieular,  we  eonsidered  the  situation 
where  a  =  1.8  em,  b  =  0.56  em,  and  H°  =  (0,  l//io)  Ampere/m.  This  is  an  ellipse 
with  an  eeeentrieity  of  0.831  and  an  aspeet  ratio  of  3.24.  We  used  a  reetangular 
simulation  domain  with  horizontal  dimension  w  =  30  em  and  vertieal  dimension 
h  =  15  em.  Eor  the  ellipse,  a  eontrived  material  was  used,  with  a  relative  perme¬ 
ability  yU  =  3.0  and  a  eonstant  isotropie  eleetrieal  eonduetivity  of  a  =  10^  S/m.  Sinee 
all  forees  and  heating  are  negleeted  in  transient  magneties,  the  density  and  tempera¬ 
ture  of  the  material  are  fixed  at  their  arbitrary  initial  values.  Eor  these  eonditions,  the 
exaet  solution  given  in  Equation  17  is  prediets  an  equilibrium  magnetie  induetion 
Bi  =  0  and  B2  =  1.18657  T  in  the  ellipse  interior. 
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Fig.  2:  Zoomed-in  view  of  meshes  used  for  ALEGRA  simulations:  body-fitted  (left) 
and  regular  (right). 


6.1  Spatial  Discretization 

Two  mesh  configurations  were  used:  (1)  a  body-fitted  irregular  mesh  representing 
an  ideal  best-case  simulation  with  no  multimaterial  elements,  and  (2)  a  regular  rect¬ 
angular  mesh  with  multimaterial  elements  representing  a  more  typical  production 
use  case  for  ALEGRA.  In  the  latter  configuration,  the  elliptical  shape  is  captured 
only  in  volume  fraction  information.  The  2  mesh  configurations  are  shown  in  Eig. 
2,  with  the  view  zoomed  in  to  show  detail  in  the  region  of  the  ellipse.  Mesh  biasing 
was  used  in  the  regular  case  to  coarsen  the  mesh  at  large  distances  exterior  to  the 
ellipse  and  concentrate  computational  effort  within  it  and  nearby;  the  regular  mesh 
bias  aspect  ratio  was  limited  to  a  maximum  of  4. 

The  meshes  shown  in  Eig.  2  were  chosen  to  provide  roughly  the  same  number  of  el¬ 
ements  inside  the  ellipse.  This  was  achieved  by  assigning  N  =  32  elements  around 
the  perimeter  of  the  ellipse  in  the  body-fitted  case,  and  =  40  elements  across  the 
major  (horizontal)  axis  in  the  regular  case.  However,  for  the  regular  mesh,  roughly 
10  times  more  elements  were  required  to  discretize  the  remainder  of  the  domain 
than  for  the  body-fitted  mesh,  even  with  mesh  biasing.  This  made  the  regular-mesh 
calculations  much  more  costly  than  the  body-fitted  mesh  calculations. 

6.2  Time  Discretization 

The  “transient  magnetics”  capability  in  ALEGRA  uses  an  implicit  solver,  and  there¬ 
fore  does  not  place  stringent  restrictions  on  the  choice  of  timestep  size.  Eor  this 
study,  the  time  evolution  of  the  system  during  relaxation  to  equilibrium  is  of  inter¬ 
est,  so  a  very  short,  fixed  timestep  dt  =  2  fis  was  still  used. 

The  time  required  for  the  system  to  reach  equilibrium  can  be  approximated  using 
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a  scaling  relationship  that  arises  from  the  1 -dimensional  (ID)  magnetic  diffusion 
equation,^  giving  a  magnetic  diffusion  time  ^  where  f  is  a  characteristic 
length  scale.  Using  the  shortest  dimension  of  the  ellipse  (b  =  0.56  cm),  we  obtain 
Trf  ~  1.2  ms  as  the  time  for  one  e-fold  increase  of  the  field  soaking  into  the  ellipse. 
Therefore,  the  simulations  were  run  out  to  a  termination  time  of  0.01  seconds  (s), 
to  capture  the  full  relaxation  to  equilibrium  by  observing  8  e-foldings. 

6.3  Solver 

The  implicit  iterative  solver  used  for  transient  magnetics  in  ALEGRA  is  the  AztecOO 
solver^^  from  the  Trilinos  package.'^  The  algebraic  multigrid  (“multilevel”)  capabil¬ 
ities  included  in  the  solver  were  enabled  here,  so  that  maximum  throughput  could 
be  obtained  on  the  larger  calculations,  and  the  actual  conditions  used  in  typical 
production  MHD  calculations  with  ALEGRA  could  be  reproduced.  Through  the 
AztecOO  interface,  a  conjugate  gradient  solver  type  was  specified,  with  an  iteration 
stopping  tolerance  of  10“^®. 

6.4  Initial  and  Boundary  Conditions 

The  ALEGRA  simulations  were  initialized  to  the  configuration  shown  in  Eig.  1, 
with  zero  magnetic  field  (and  thus  zero  current)  everywhere  in  the  domain.  The 
system  becomes  magnetized  because  of  a  uniform  y-oriented  source  magnetic  field 
imposed  on  the  left  and  right  boundaries,  which  is  assumed  to  be  switched  on  in¬ 
stantaneously  at  time  zero.  The  component  of  the  magnetic  field  if*  tangent  to  the 
mesh  is  fixed  on  the  left  and  right  boundaries,  so  that  over  the  course  of  the  simu¬ 
lation,  the  solution  relaxes  to  an  equilibrium  state,  which  is  equivalent  to  the  equi¬ 
librium  condition  described  in  Section  2.  On  the  top  and  bottom  boundaries,  the 
tangent  component  of  if*  is  constrained  to  be  zero,  so  that  field  lines  are  normal  to 
the  boundary  there. 

Explicit  boundary  conditions  on  the  material  interfaces  are  not  needed.  The  mate¬ 
rial  discontinuity  is  represented  in  the  calculation  as  a  discontinuity  in  a  and  /i.  In 
multimaterial  elements,  mean  values  of  these  properties  are  used. 
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Fig.  3:  Configuration  of  field  lines  before  a)  and  after  h)  after  the  equilibrium 
state  is  reached,  for  the  body-fitted  (left)  and  regular  (right)  mesh  types.  Simulation 
times:  a)  0.2  ms,  h)  3.0  ms. 

7.  Time  Evolution  of  the  System 

In  the  ALEGRA  simulations,  the  evolving  variables  are  the  in-plane  magnetic  in¬ 
duction  vector  B\  the  in-plane  magnetic  field  H\  the  out-of-plane  magnetic  vector 
potential  A,  and  the  out-of-plane  electric  current  density  J.  The  latter  2  vectors 
are  constrained  to  have  only  a  2;  component,  so  are  treated  by  the  code  as  scalar 
quantities. 

The  code  also  computes  a  total  electric  current  J,  which  is  a  positive  quantity  ob¬ 
tained  from  surface  integrals  of  J.  The  total  surface  integral  is  zero  by  construction; 
this  can  be  seen  by  taking  the  divergence  of  Equation  42.  The  magnitude  of  opposed 
currents  running  transiently  in  the  conductor  normal  to  the  plane  is  obtained  by  iso¬ 
lating  the  positive  and  negative  contributions  to  the  total  surface  integral,  which 
should  be  equal. 

Over  the  course  of  the  simulation,  the  current  /  and  time  derivatives  of  H\  B\ 
and  A  all  relax  toward  zero  as  magnetic  diffusion  proceeds  toward  a  “saturated” 
or  equilibrium  field  configuration.  Magnetic  field  lines  computed  by  ALEGRA  be¬ 
fore  and  after  equilibrium  is  reached  are  shown  for  both  mesh  types  in  Eig.  3.  The 
lines  are  obtained  conveniently  in  practice  by  plotting  contours  of  A.  In  the  system 


16 


considered  here,  with  /r  7^  /rO,  the  equilibrium  field  configuration  is  nonuniform. 
As  predicted  in  Section  2,  it  is  only  uniform  in  the  interior  of  the  inclusion  late  in 
time.  The  uniform  interior  field  configuration  has  visibly  not  been  reached  yet  by 
t  =  0.6  /iS,  as  seen  in  Fig.  3a. 

The  evolution  of  the  magnitude  of  5*  and  J  in  the  ALEGRA  simulations  is  de¬ 
picted  in  Figs.  4  and  5.  We  observe  that  the  field  “soaks”  into  the  ellipse  nonuni- 
formly.  Before  the  equilibrium  value  is  reached,  the  local  magnitude  of  the  field 
first  transiently  rises  to  much  larger  values  on  the  extreme  horizontal  limbs  of  the 
ellipse.  The  field  magnitude  near  the  center  of  the  ellipse  remains  small  until  these 
transients  decay,  leading  to  a  uniform  vertically  oriented  interior  field. 

Fig.  4  shows  results  for  the  body-fit  mesh,  and  Fig.  5  shows  results  for  the  regular 
mesh.  Mesh  lines  are  superimposed  in  both  cases  in  the  region  of  the  ellipse.  In 
the  regular-mesh  case,  mesh  lines  are  plotted  for  material  volume  fraction  /  >  0.5, 
with  the  estimated  contour  of  /  =  0.5  appearing  as  well.  It  should  be  noted  that  this 
contour  does  not  represent  a  mesh  line,  but  a  post-processing  contour  estimate. 

Comparing  the  body-fitted  and  regular-mesh  results,  we  see  that  there  are  advan¬ 
tages  and  disadvantages  to  both  approaches.  A  much  stronger  discontinuity  in  the 
solution  variables  is  possible  at  the  material  interface  in  the  body-fitted  case,  allow¬ 
ing  greater  accuracy.  This  greater  accuracy  is  enormous  for  this  coarse  mesh  size. 
However,  the  elements  are  much  more  distorted,  allowing  other  errors  to  enter. 

The  primary  advantage  of  the  body-fit  approach  is  that  a  much  smaller  total  num¬ 
ber  of  elements  is  required  for  this  simulation  domain:  only  about  900,  compared 
to  over  11,000  for  the  regular  case.  This  allows  computational  effort  to  be  concen¬ 
trated  in  the  region  of  interest  (the  ellipse).  The  smaller  element  count  results  in 
a  much  shorter  simulation  turnaround  time:  the  body-fitted  case  runs  to  comple¬ 
tion  (in  serial,  for  this  coarse  mesh  size)  in  about  15%  of  the  time  required  for  the 
regular-mesh  case  (see  table  in  Section  8). 

The  accuracy  per-unit-compute-time  is  enormously  greater  with  a  body-fitted  mesh. 
However,  in  most  production  calculations,  a  body-fitted  mesh  is  not  a  viable  option. 
Therefore,  we  have  retained  both  mesh  types  in  this  study. 
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0.2  ms 


0.6  ms 


1.4  ms 


3  ms 


Fig.  4:  Simulation  evolution  for  body-fitted  mesh:  magnitude  (left)  and  J  mag¬ 

nitude  (right),  with  mesh  lines  in  the  ellipse  superimposed.  Simulation  times:  a)  0.2 
ms,  h)  0.6  ms,  c)  1.4  ms,  d)  3  ms.  Blue-red  eolorbar  for  spans  0  <  |i?*|  <  2  T. 
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(d) 
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Fig.  5:  Simulation  evolution  for  regular  mesh:  S*  magnitude  (left)  and  J  magnitude 
(right),  with  mesh  lines  in  the  ellipse  superimposed.  Simulation  times:  a)  0.2  ms,  h) 
0.6  ms,  c)  1.4  ms,  d)  3  ms.  Blue-red  colorbar  for  spans  0  <  \B^\  <  2  T. 
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Fig.  6:  Spatial  distribution  of  deviation  from  exact  solution  for  B2  at  equilibrium 
for  body-fitted  mesh  with  iV  =  40. 

8.  Solution  Verification 

The  interior  solution  at  f  =  0.01  s  was  compared  to  the  equilibrium  solution  given 
by  Equation  17,  which  predicts  B2  =  1.18657  T.  In  the  case  of  the  regular  mesh, 
only  elements  with  /  >  0.5  were  considered  in  the  comparison.  Since  B2  is  a  de¬ 
rived  variable  in  ALEGRA,  and  numerical  differentiation  is  involved  in  computing 
its  value,  some  additional  inaccuracy  is  involved  in  the  comparison.  The  fundamen¬ 
tal  variable  in  these  calculations  is  A,  but  comparison  was  not  done  for  A  here, 
because  in  production  calculations,  the  vector  B^  is  typically  of  principal  interest. 
Nevertheless,  even  for  B2,  we  expect  a  reasonable  qualitative  match  to  the  exact 
solution,  and  at  least  first-order  convergence  toward  the  solution  under  spatial  mesh 
refinement,  based  on  the  methods  implemented  in  AEEGRA,  and  the  mesh  types 
used  here.'’ 

8.1  Coarse-Mesh  Comparison 

Comparison  of  the  equilibrium  distribution  of  B2  in  the  ellipse  for  the  2  mesh  types 
shows  that  for  the  =  40  case,  which  is  a  very  coarse  mesh  requiring  only  minutes 
of  serial  compute  time,  AEEGRA  matches  the  exact  solution  from  Equation  17  to 
within  fractions  of  a  percent.  The  distributions  computed  on  the  =  40  regular 
and  body-fitted  mesh  types  are  shown  with  the  exact  solution  in  Eigs.  6  and  7.  The 
field  is  not  perfectly  uniform,  but  deviates  from  the  uniform  exact  solution  by  only 
0.02%  for  the  body-fitted  case  and  by  only  0.1%  for  the  regular  case. 
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Fig.  7:  Spatial  distribution  of  deviation  from  exact  solution  for  B2  at  equilibrium 
for  regular  mesh  with  =  40. 
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The  largest  deviations  appear  on  the  horizontal  limbs  of  the  ellipse,  particularly  in 
the  regular-mesh  case.  This  suggests  that  greater  inaccuracies  are  likely  for  elliptical 
inclusions  with  greater  aspect  ratios  than  3.24  (see  Section  8  below).  The  deviations 
are  also  much  larger  for  the  regular  mesh  than  the  body-fitted  mesh  at  this  mesh  size 
because  of  the  effect  of  mixed  elements  near  the  ellipse  surface.  At  the  center  of 
the  ellipse,  the  field  magnitude  exceeds  the  analytic  value.  Further  analysis  shows 
that  the  mean  value  of  the  computed  field  in  the  ellipse  does  slightly  exceed  the 
analytical  value,  as  suggested  by  Figs.  6  and  7. 

The  global  deviation  from  the  exact  solution  was  computed  by  summing  the  square 
of  the  local  deviation  over  all  locations  inside  the  ellipse  (/  >  0.5),  and  normalizing 
by  the  tally  of  elements  inside  the  ellipse.  Taking  the  root  then  yields  a  root-mean- 
square  (RMS)  measure  of  the  error  in  the  computed  B  —  2  distribution.  Normalizing 
again  by  the  value  of  the  exact  solution  yields  a  fractional  error. 

This  fractional  error  metric  is  plotted  against  time  in  Fig.  8  for  the  =  40  case. 
We  observe  again  that  the  body-fitted  result  is  significantly  more  accurate  at  this 
mesh  size.  Further,  we  see  that  the  error  trend  in  both  cases  goes  through  a  min¬ 
imum  during  the  transient  phase,  marking  the  moment  when  the  average  value  of 
the  field  in  the  initially  ellipse  first  surpasses  the  analytic  value.  The  error  reaches  a 
plateau  after  about  f  =  5  ms.  Also  at  about  f  =  5  ms,  the  current  ceases  to  change 
significantly  in  time,  reaching  its  final  small  but  nonzero  value. 
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Fig.  8:  Time  evolution  of  RMS  deviation  of  ALEGRA  solution  relative  to  exaet 
solution  for  equilibrium  eondition  at  =  40. 


8.2  Convergence  Under  Spatial  Refinement 

Spatial-refinement  eonvergenee  studies  were  also  performed  for  both  mesh  types, 
starting  from  the  eoarse  mesh  shown  in  the  foregoing  diseussion.  Previous  veri- 
fieation  work  for  transient  magnetics  in  ALEGRA  has  shown  good  convergence 
behavior  in  3-D.^^’^’  Here,  the  2-D  mesh  interval  was  halved  5  times,  yielding  6 
resolution  levels  at  =  40,  80, 160,  320,  640,  and  1,  280  elements  per  perimeter 
(body-fitted  case)  or  elements  per  major  axis  (regular  case)  for  comparison  to  the 
exact  solution,  using  the  fractional  RMS  error  metric  just  described. 

The  error  computed  in  this  way  from  the  series  of  simulations  was  observed  to 
decay  monotonically  at  approximately  first  order  relative  to  the  mesh  interval.  The 
error  convergence  trends  are  shown  in  Fig.  9.  The  trends  are  plotted  against  the 
total  number  of  elements  in  the  ellipse  for  each  simulation,  but  convergence  rates 
are  computed  against  the  value  of  N.  We  find  average  convergence  rates  of  0.9  for 
the  body-fitted  meshes  and  1.3  for  the  regular  meshes. 

We  see  that  the  transient  magnetics  module  in  ALEGRA  computes  solutions  that 
converge  smoothly  toward  the  analytic  solution  at  the  expected  rate.  ALEGRA 
computes  magnetic  diffusion  into  an  ellipse  with  linear  magnetic  permeability  with 
demonstrated  accuracy,  and  the  expected  convergence  behavior. 
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Fig.  9:  Results  of  spatial  mesh  refinement  eonvergenee  study  for  an  elliptieal  inelu- 
sion,  showing  fractional  error  relative  to  exact  solution. 


Under  spatial  refinement,  it  becomes  clear  that  the  regular  mesh  approach  with 
mixed  elements  is  not  as  disadvantaged  relative  to  body-fitted  meshes  as  the  =  40 
results  indicate.  Although  the  error  magnitudes  are  larger  for  the  regular  mesh  with 
mixed  elements,  the  solution  computed  on  these  meshes  converges  at  a  much  higher 
rate.  Its  convergence  rate  also  improves  as  the  mesh  is  refined.  At  the  very  largest 
meshes,  the  error  magnitudes  are  comparable.  Thus,  there  is  a  tradeoff  between 
regular  and  body-fitted  meshes,  exchanging  accuracy  for  convergence. 

8.3  Aspect  Ratio  Effect 

Since  the  performance  of  ALEGRA  for  a  range  of  aspect  ratios  is  also  of  interest, 
the  spatial  mesh  refinement  study  was  repeated  for  2  more  geometries:  a  circular 
inclusion,  and  an  elliptical  inclusion  with  an  aspect  ratio  of  10.  In  the  case  of  a  cir¬ 
cular  inclusion,  the  dimensions  a  =  b  =  1  cm  were  used,  with  domain  dimensions 
w  =  h  =  20  cm.  Otherwise  the  parameterization  and  analysis  were  the  same.  The 
exact  solution  for  this  case  predicts  an  equilibrium  vertical  magnetic  induction  of 
i?2  =  1-5  T  in  the  interior  of  the  circle.  Body-fitted  and  regular  meshes  were  used, 
with  the  same  series  of  mesh  refinements. 
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Fig.  10:  Results  of  spatial  mesh  refinement  eonvergenee  study  for  a  circular  inelu¬ 
sion,  showing  fraetional  error  relative  to  exaet  solution. 


Results  for  the  eireular  ease  are  shown  in  Fig.  10.  Not  surprisingly,  the  error  mag¬ 
nitudes  are  for  the  most  part  lower  for  this  simpler  ease,  at  eaeh  value  of  N.  We 
also  again  observe  larger  error  magnitudes  but  greater  eonvergenee  rates  for  the 
regular  mesh  than  for  the  body-fitted  mesh.  Further,  the  disparity  in  eonvergenee 
rates  also  remains.  The  regular-mesh  ease  eonverges  at  a  mean  rate  of  1 .5,  while  the 
body-fitted  ease  at  only  0.7. 

To  simulate  an  ellipse  with  aspeet  ratio  AR  =  10,  dimensions  a  =  3.16  em  and 
6  =  0.316  em  were  used,  with  domain  dimensions  w  =  50  em  and  h  =  20  em.  This 
eorresponds  to  an  eeeentrieity  of  0.949.  The  exaet  solution  for  this  ease  prediets 
an  equilibrium  vertieal  magnetie  induetion  of  B2  =  1.06452  T  in  the  interior  of 
the  ellipse.  For  this  large  aspeet  ratio,  the  meshing  teehnique  used  to  produee  the 
elliptieal  body-fitted  meshes  was  not  operable  due  to  the  large  mesh  deformations 
required.  Henee,  no  body-fit-mesh  eases  were  studied.  Further,  in  the  regular-mesh 
ease,  a  lower- aeeuraey  geometry  insertion  method  (“diatom,  insert  UDS”  rather 
than  “insert  UDS-i-”)  had  to  be  used,  beeause  of  the  high  eurvature  of  the  shape. 

eonvergenee  results  for  the  AR  =  10  ease  are  shown  in  Fig.  11,  along  with  the  re¬ 
sults  from  the  AR  =  3.24  and  AR  =  1  eases  diseussed  above  with  regular  meshes. 
Here  the  fraetional  error  is  plotted  as  a  funetion  of  N,  not  the  number  of  elements 
in  the  ellipse. 
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Fig.  11:  Results  of  spatial  mesh  refinement  eonvergenee  study  for  elliptieal  inelu- 
sions  with  various  aspeet  ratios,  showing  fraetional  error  relative  to  the  respeetive 
exact  solution  as  a  function  of  N,  for  regular-mesh  simulations. 


For  AR  =  10,  the  computed  solution  converges  less  smoothly,  and  the  error  mag¬ 
nitudes  are  larger  than  for  smaller  aspect  ratios.  However,  the  solution  converges  at 
a  mean  rate  of  1.2,  which  still  compares  favorably  with  the  convergence  rates  for 
the  circular  inclusion  (1.5)  and  the  original  elliptical  inclusion  (1.3).  This  suggests 
that  with  the  regular-mesh  setup,  ALEGRA’s  accuracy  for  this  problem  should  be 
preserved  even  for  ellipses  with  very  great  eccentricity. 

8.4  Compute  Times 

As  a  final  note,  some  diagnostic  data  from  the  simulations  for  aspect  ratio  3.24 
are  provided  in  the  table,  to  illustrate  the  disparity  in  computational  cost  for  the 
body-fitted  and  regular  mesh  scenarios.  This  includes  the  number  of  elements  in  the 
mesh,  the  number  of  compute  cores  used  for  each  simulation,  and  the  “wall-clock” 
time  required  for  the  ALEGRA  calculation  to  run.  The  regular-mesh  simulations 
required  10  x  more  elements  at  each  value  of  N,  resulting  in  a  much  lower  rate  of 
throughput,  due  to  higher  AztecOO  solve  times. 

Despite  the  apparent  lower  rate  of  throughput  for  the  regular-mesh  cases,  there  are  2 
mitigating  factors.  The  body-fitted  cases  all  carried  fewer  mesh  elements  per  com¬ 
pute  core,  by  about  a  factor  of  2.  Also,  the  body-fitted  cases  all  required  additional 
compute  time  to  generate  the  mesh.  The  mesh-generation  process  does  not  lend 
itself  to  parallelization,  and  for  the  larger  meshes  required  more  time  than  the  ac- 
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N 

^elem 

Body-fit 

^elem 

Regular 

^cores 

Body-fit 

^cores 

Regular 

Wall  time  (s) 
Body-fit 

Wall  time  (s) 
Regular 

40 

1.30x10^ 

1.12x10^ 

2 

1 

9.20x10^ 

8.30x102 

80 

4.80x10^ 

4.58x10^ 

2 

4 

2.26x10^ 

1.29x10® 

160 

1.84x10^ 

1.85x10® 

4 

8 

5.53x102 

3.76x10® 

320 

7.20x10^ 

7.42x10® 

8 

16 

1.27x10® 

9.32x10® 

640 

2.85  xlO^ 

2.98x10® 

16 

32 

3.61x10® 

2.22x10^ 

1280 

1.13x10® 

1.19x10^ 

32 

64 

9.31x10® 

6.62x10^ 

Table:  Diagnostic  data  from  simulations  run  in  spatial-refinement  convergence 
study. 


tual  ALEGRA  calculation.  The  expense  of  this  operation  mainly  has  to  do  with  the 
Laplacian  smoothing  operation  that  had  to  be  used  to  improve  the  symmetry  of  the 
quadrilateral  mesh;  this  operation  might  not  be  required  for  triangular  meshes.  This 
is  a  significant  drawback  to  using  the  body-fitted  mesh  approach. 

9.  Conclusions 

This  study  has  yielded  several  outcomes  that  should  be  useful  for  those  interested  in 
understanding  the  quasi-static  behavior  of  magnetic  fields  in  a  medium  that  includes 
conducting,  permeable  materials. 

First,  a  flexible  and  convenient  analytic  solution  has  been  derived  that  describes  the 
equilibrium  field  in  the  interior  of  an  inclusion  in  such  a  medium  that  has  a  generic 
shape:  an  ellipse.  Although  the  solution  is  rather  complex  outside  the  ellipse,  inside 
the  ellipse,  it  has  a  simple,  closed  form  which  lends  itself  well  to  verification  of 
computed  solutions  by  both  visual  inspection  and  convergence  analysis.  A  similar 
static  solution  is  valid  in  the  electrostatic  context. 

Second,  it  has  been  shown  that  the  equilibrium  state  can  be  reached  by  transient 
means  via  computation  with  ALEGRA,  which  allows  the  magnetic  field  to  evolve 
by  solving  an  initial-boundary- value  problem  that  includes  linear  magnetic  perme¬ 
ability  in  the  inclusion.  The  time  to  reach  the  equilibrium  state  matches  the  mag¬ 
netic  diffusion  time  expected  for  the  system  based  on  its  geometry  and  the  material 
properties.  This  approach  opens  a  novel  venue  for  solving  classical  problems  of 
magnetostatics  and  electrostatics  of  dielectrics. 

Third,  it  has  been  demonstrated  that  the  equilibrium  state  computed  by  ALEGRA 
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for  this  system  reliably  reproduces  the  analytic  solution  quite  well,  with  errors  on 
the  order  of  10“^  to  10“^  for  reasonable  compute  times.  Further,  the  solution  con¬ 
verges  under  spatial  refinement,  and  the  convergence  properties  are  preserved  even 
for  ellipse  aspect  ratios  up  to  10. 

Finally,  the  results  show  that  it  is  not  necessarily  more  advantageous  to  use  body- 
fitted  meshes  when  simulating  this  system  in  ALEGRA  on  a  quadrilateral  mesh. 
For  coarse  meshes,  the  body-fitted  approach  does  provide  much  greater  accuracy 
and  faster  throughput.  But  as  the  mesh  is  refined,  the  advantage  in  accuracy  slowly 
vanishes,  and  the  cost  of  generating  the  body-fitted  mesh  becomes  excessive,  and 
even  prohibitive  for  large  ellipse  aspect  ratios.  At  all  mesh  sizes,  the  convergence 
rate  on  regular  meshes  surpasses  that  of  body-fitted  meshes. 

This  work  lays  the  foundation  for  future  extension  into  other  areas  of  interest.  An 
analytic  solution  in  three  dimensions  for  an  ellipsoidal  inclusion  is  also  possible, 
which  can  be  compared  to  computed  solutions  from  ALEGRA  to  test  its  3-D  tran¬ 
sient  magnetics  capabilities.  Electromagnetic  forces,  motion,  and  thermodynamic 
effects  could  also  be  incorporated,  without  introducing  significant  difficulties  in 
ALEGRA.  All  of  these  are  areas  that  could  be  explored  in  future  analysis  with 
fruitful  results. 
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List  of  Symbols 


a  -  Horizontal  ellipse  semi-major  axis 
AR  -  Ellipse  aspeet  ratio 
Bi  -  Magnetie  induetion  veetor 
b  -  Vertical  ellipse  semi-minor  axis 
c  -  Speed  of  light 
Ei  -  Electric  field  vector 
Hi  -  Magnetic  field  vector 
-  Magnetic  field  at  infinity 
h  -  Vertical  dimension  of  simulation  domain 
J*  -  Current  density  vector 
yU  -  Material  magnetic  permeability 
/io  -  Magnetic  permeability  of  free  space 
Ui  -  Unit  normal  vector 

N  -  Number  of  elements  spanning  ellipse  major  axis  or  perimeter 
a  -  Isotropic  material  electrical  conductivity  scalar 
cr*-^  -  Electrical  conductivity  tensor 
t  -  time 

Td  -  Magnetic  diffusion  time 
Tj  -  Unit  tangent  vector 
0  -  Eogarithmic  potential 
w  -  Horizontal  dimension  of  simulation  domain 
Zijk  -  Eevi-Civita  tensor 
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